This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for Pacuta larvae.
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car')
if ("lme4" %in% rownames(installed.packages()) == 'FALSE') install.packages('lme4')
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest')
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales')
if ("cowplot" %in% rownames(installed.packages()) == 'FALSE') install.packages('cowplot')
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2')
if ("effects" %in% rownames(installed.packages()) == 'FALSE') install.packages('effects')
if ("emmeans" %in% rownames(installed.packages()) == 'FALSE') install.packages('emmeans')
if ("multcomp" %in% rownames(installed.packages()) == 'FALSE') install.packages('multcomp')
#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('lme4')
library('lmerTest')
library('scales')
library('cowplot')
library('effects')
library('emmeans')
library('multcomp')
Load data from LoLinR.
PRdata<-read.csv("Pacu2021/Output/Respiration/PR/oxygen_P_R_calc.csv") #load data
Format data columns.
#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]
#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)
Calculate a P:R ratio using gross photosynthesis and respiration calculated as GP:R.
PRdata$ratio<-abs(PRdata$GP.nmol.org.min)/abs(PRdata$R.nmol.org.min) #calculate ratio with absolute values
View outliers. Remove outliers for P:R data.
boxplot(PRdata$R.nmol.org.min)
boxplot(PRdata$P.nmol.org.min)
boxplot(PRdata$GP.nmol.org.min)
boxplot(PRdata$ratio)
PRdata<-PRdata%>%filter(ratio < 4) #filter out outlier for PR ratio
boxplot(PRdata$ratio)
Plot data as a box plot.
r_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
PRdata$Date<-as.factor(PRdata$Date)
r_plot1<-PRdata %>%
group_by(Temperature, as.factor(Date))%>%
ggplot(., aes(x = as.factor(Temperature), y = R.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.2, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot1
Plot data as a mean dot plot.
r_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(R.nmol.org.min), sd=sd(R.nmol.org.min), N=length(R.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.1, 0.06), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot_b
For photosynthesis filter out runs 7 and 8, these runs were outliers that didn’t work due to equipment issues.
Plot data as a box plot.
p_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot
p_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = P.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot1
Plot data as a mean dot plot.
p_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(P.nmol.org.min), sd=sd(P.nmol.org.min), N=length(P.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot_b
Plot data as a box plot.
gp_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot
gp_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = GP.nmol.org.min, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.02, 0.33), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot1
Plot data as a mean dot plot.
gp_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(GP.nmol.org.min), sd=sd(GP.nmol.org.min), N=length(GP.nmol.org.min), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.005, 0.15), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot_b
Plot data as a box plot.
pr_plot<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot
pr_plot1<-PRdata %>%
group_by(Temperature)%>%
ggplot(., aes(x = as.factor(Temperature), y = ratio, group=Temperature)) +
facet_wrap(~Date)+
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(fill=Temperature, group=Temperature), alpha=0.5)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=4, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 4), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot1
Plot data as a mean dot plot.
pr_plot_b<-PRdata %>%
group_by(Temperature)%>%
summarise(mean=mean(ratio), sd=sd(ratio), N=length(ratio), se=sd/sqrt(N))%>%
ggplot(., aes(x = as.factor(Temperature), y = mean, group=Temperature)) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_point(aes(fill=Temperature, group=Temperature), pch = 21, size=6, alpha=0.8, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se, group=Temperature), width=0, linetype="solid", position=position_dodge(0.3), size=1, color="black")+
xlab("Temperature") +
scale_fill_manual(name="Temperature", values=c("blue","orange", "red"))+
scale_color_manual(name="Temperature", values=c("blue","orange", "red"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(0, 2), labels = scales::number_format(accuracy = 0.1, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); pr_plot_b
Box plot assembled figure.
full_fig<-plot_grid(r_plot, p_plot, gp_plot, pr_plot, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")
ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_boxplot.png", plot=full_fig, dpi=500, width=18, height=6, units="in")
Mean plot assembled figure.
full_fig_b<-plot_grid(r_plot_b, p_plot_b, gp_plot_b, pr_plot_b, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")
ggsave(filename="Pacu2021/Figures/Respiration/PR/respirometry_fig_mean.png", plot=full_fig_b, dpi=500, width=18, height=6, units="in")
Change data columns to factors.
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Date<-as.factor(PRdata$Date)
PRdata$Run<-as.factor(PRdata$Run)
View respiration data.
hist(PRdata$R.nmol.org.min)
Build model and examine for respiration across temperatures.
Rmodel1<-aov(R.nmol.org.min~Temperature, data=PRdata)
summary(Rmodel1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.00617 0.003087 2.537 0.082 .
## Residuals 180 0.21903 0.001217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check for effects of day.
Rmodel1b<-aov(R.nmol.org.min~Temperature*Date, data=PRdata)
summary(Rmodel1b)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.00617 0.003087 2.562 0.080 .
## Date 2 0.00454 0.002269 1.884 0.155
## Temperature:Date 1 0.00123 0.001228 1.019 0.314
## Residuals 177 0.21326 0.001205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance.
qqPlot(residuals(Rmodel1))
## [1] 131 123
leveneTest(residuals(Rmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.0868 0.1271
## 180
Conduct post hoc test.
emm<-emmeans(Rmodel1, ~Temperature)
pairs(emm)
## contrast estimate SE df t.ratio p.value
## 27 - 31 0.01214 0.00590 180 2.056 0.1022
## 27 - 34 -0.00127 0.00688 180 -0.184 0.9815
## 31 - 34 -0.01340 0.00743 180 -1.805 0.1708
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Build model and examine for Photosynthesis across temperatures.
Pmodel1<-aov(P.nmol.org.min~Temperature, data=PRdata)
summary(Pmodel1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0157 0.007851 7.986 0.000475 ***
## Residuals 180 0.1769 0.000983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check for effect of day.
Pmodel1b<-aov(P.nmol.org.min~Temperature*Date, data=PRdata)
summary(Pmodel1b)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.01570 0.007851 8.328 0.000349 ***
## Date 2 0.01008 0.005040 5.347 0.005565 **
## Temperature:Date 1 0.00001 0.000009 0.009 0.922918
## Residuals 177 0.16685 0.000943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance.
qqPlot(residuals(Pmodel1))
## [1] 88 58
leveneTest(residuals(Pmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 7.2644 0.0009246 ***
## 180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conduct post hoc test.
emm<-emmeans(Pmodel1, ~Temperature)
pairs(emm)
## contrast estimate SE df t.ratio p.value
## 27 - 31 0.00579 0.00531 180 1.090 0.5215
## 27 - 34 0.02466 0.00618 180 3.989 0.0003
## 31 - 34 0.01888 0.00667 180 2.828 0.0144
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Build model and examine for gross photosynthesis across temperatures.
GPmodel1<-aov(GP.nmol.org.min~Temperature, data=PRdata)
summary(GPmodel1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0247 0.012359 4.241 0.0159 *
## Residuals 180 0.5246 0.002914
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check for effect of day.
GPmodel1b<-aov(GP.nmol.org.min~Temperature*Date, data=PRdata)
summary(GPmodel1b)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 0.0247 0.012359 4.395 0.0137 *
## Date 2 0.0253 0.012660 4.501 0.0124 *
## Temperature:Date 1 0.0014 0.001445 0.514 0.4745
## Residuals 177 0.4978 0.002812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance.
qqPlot(residuals(GPmodel1))
## [1] 41 131
leveneTest(residuals(GPmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2835 0.2796
## 180
Conduct post hoc test.
emm<-emmeans(GPmodel1, ~Temperature)
pairs(emm)
## contrast estimate SE df t.ratio p.value
## 27 - 31 -0.00635 0.00914 180 -0.695 0.7666
## 27 - 34 0.02593 0.01065 180 2.436 0.0417
## 31 - 34 0.03228 0.01149 180 2.809 0.0152
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Build model and examine for PR ratio across temperatures.
PRmodel1<-aov(ratio~Temperature, data=PRdata)
summary(PRmodel1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 2.89 1.4471 4.087 0.0184 *
## Residuals 180 63.73 0.3541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check for effect of day.
PRmodel1b<-aov(ratio~Temperature*Date, data=PRdata)
summary(PRmodel1b)
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 2 2.89 1.4471 4.339 0.01447 *
## Date 2 4.62 2.3078 6.920 0.00128 **
## Temperature:Date 1 0.08 0.0837 0.251 0.61701
## Residuals 177 59.03 0.3335
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance.
qqPlot(residuals(PRmodel1))
## [1] 159 94
leveneTest(residuals(PRmodel1)~Temperature, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 4.1833 0.01676 *
## 180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conduct post hoc test.
emm<-emmeans(PRmodel1, ~Temperature)
pairs(emm)
## contrast estimate SE df t.ratio p.value
## 27 - 31 0.189 0.101 180 1.873 0.1495
## 27 - 34 0.312 0.117 180 2.660 0.0230
## 31 - 34 0.124 0.127 180 0.975 0.5936
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Generate summary of all respiration data.
summary<-PRdata%>%
group_by(Temperature)%>%
filter(R.nmol.org.min<0)%>%
filter(P.nmol.org.min>0)%>%
filter(GP.nmol.org.min>0)%>%
summarise(N=length(R.nmol.org.min),
Mean_R=mean(R.nmol.org.min),
SE_R=sd(R.nmol.org.min)/sqrt(length(R.nmol.org.min)),
Mean_P=mean(P.nmol.org.min),
SE_P=sd(P.nmol.org.min)/sqrt(length(P.nmol.org.min)),
Mean_GP=mean(GP.nmol.org.min),
SE_GP=sd(GP.nmol.org.min)/sqrt(length(GP.nmol.org.min)),
Mean_PR=mean(ratio),
SE_PR=sd(ratio)/sqrt(length(ratio)))
summary
## # A tibble: 3 × 10
## Temperature N Mean_R SE_R Mean_P SE_P Mean_GP SE_GP Mean_PR
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 27 81 -0.0728 0.00345 0.0577 0.00358 0.130 0.00531 1.91
## 2 31 54 -0.0828 0.00374 0.0485 0.00367 0.131 0.00580 1.66
## 3 34 32 -0.0738 0.00690 0.0308 0.00295 0.105 0.00799 1.55
## # … with 1 more variable: SE_PR <dbl>
summary%>%
write_csv(., "Pacu2021/Output/Respiration/PR/mean_respiration.csv")